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Abstract. A fast multilevel algorithm based on directionally scaled 
tensor-product Gaussian kernels on structured sparse grids is proposed 
for interpolation of high-dimensional functions and for the numerical in¬ 
tegration of high-dimensional integrals. The algorithm is based on the 
recent Multilevel Sparse Kernel-based Interpolation (MLSKI) method 
(Georgoulis, Levesley & Subhan, SIAM J. Sci. Comput., 35(2), pp. A815- 
A831, 2013), with particular focus on the fast implementation of Gaussian- 
based MLSKI for interpolation and integration problems of high-dimen¬ 
sional functions / : [0, l] d — > R, with 5 < d < 10. The MLSKI interpo¬ 
lation procedure is shown to be interpolatory and a fast implementation 
is proposed. More specifically, exploiting the tensor-product nature of 
anisotropic Gaussian kernels, one-dimensional cardinal basis functions on 
a sequence of hierarchical equidistant nodes are precomputed to machine 
precision, rendering the interpolation problem into a fully parallelisable 
ensemble of linear combinations of function evaluations. A numerical 
integration algorithm is also proposed, based on interpolating the (high¬ 
dimensional) integrand. A series of numerical experiments highlights the 
applicability of the proposed algorithm for interpolation and integration 
for up to 10-dimensional problems. 


1 Introduction 

Approximation, interpolation and numerical integration of high-dimensional func¬ 
tions have numerous applications in applied sciences, ranging from mathematics, 
statistics, physics and engineering, to economics and finance. 

For low dimensional multivariate functions / : fl —> R, fl C M d , with 2 < d < 
4, classical interpolation and integration (cubature) methods are able to deliver 
accurate and efficient results to high accuracy, see e.g., 025 ] and the references 
therein. For d > 5 , standard algorithms are challenged by the, so-called, curse of 
dimensionality. Indeed, standard polynomial, spline or kernel-based interpolation 
methods of a function / : [0, l] d -> R, on a d-dimensional grid with N nodes 
in each direction, typically require N d function evaluations to achieve rates of 
convergence of order N ~ a , for some a > 0 independent of d; a manifestation of 
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the curse of dimensionality, as the interpolation problem becomes exponentially 
more computationally intensive with d. 

To address this issue, a number of methods based on either randomness, 
especially for integration (Monte-Carlo, quasi-Monte-Carlo, multilevel Monte- 
Carlo; see, e.g., I.M'ill'iloljnl and the references therein) or on Smolyak/sparse 
grid/hyperbolic cross/boolean interpolation-type constructions (see, e.g., [241812612314] 
and the references therein) have appeared in the literature during the last five 
decades. Monte-Carlo-type approaches are, generally speaking, robust and easy 
to implement and their accuracy is independent of the problem dimension d, 
at the cost of typically slow convergence in a probabilistic sense m ■ Sparse 
grids and hyperbolic crosses are applicable when the function / has sufficiently 
regular mixed derivatives; in such cases, they can deliver, typically algebraic, 
rates of convergence, independent of d with typical complexity of IV (log IV) d_1 , 
where N is the number of points in each direction. More recently, “p-version” 
sparse grids, based on polynomial interpolation on Chebyshev or Clenshaw- 
Curtis sparse tensor-product type grids have been proposed for the calculation 
of high-dimensional integrals arising in the numerical approximation of elliptic 
boundary-value problems with random coefficients mm- 

Another approach, that can be positioned within the second class of meth¬ 
ods described above is the recent Multilevel Sparse Kernel-based Interpolation 
(MLSKI) method T2] . MLSKI method is based on the concept of combination 
(also known as d-dimensional boolean interpolation [8] ) of smaller interpolation 
problems on directionally uniform nodes, whose union coincides with the sparse 
grid node systems produced by hierarchical linear splines. Each interpolation 
sub-problem is evaluated using anisotropically scaled versions of translation of 
standard kernels/radial basis functions. All interpolation problems are then lin¬ 
early combined appropriately, leading to the Sparse Kernel-based Interpolation 
method (SKI). To further exploit the approximation power of the SKI inter- 
polants, which use naturally nested sparse grid nodes, the interpolation proce¬ 
dure is performed in a multilevel fashion, starting from a coarse SKI interpolant 
so of the function /, and then continuing by interpolating the residual / — So on 
the next level sparse grid, and so on. We stress that, despite the fact that the 
MLSKI method is applied on sparse grids arising from uniform full grids, it has 
not been observed to be susceptible to classical Runge-type instabilities. Hence, 
MLSKI can be viewed as an alternative to “p-version” sparse grids, by offering 
comparable convergence rates on “uniform” (and, therefore, hierarchical) sparse 
grids. 

This work is concerned with further investigating the MLSKI approach both 
theoretically and numerically. More specifically, focusing on MLSKI using Gaus¬ 
sian kernels, we show that at each level of SKI, the method indeed, interpolates 
the function /. To achieve numerical exactness on all the nodes, the MLSKI 
approach it is sufficient to use tensor-product-type kernels, i.e., Gaussians. This 
also leads naturally to the fast implementation of the sub-problem computations 
of MLSKI 14], each of which can be computed completely independently. In 
particular, exploiting the tensor-product nature of anisotropic Gaussian kernels, 
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one-dimensional cardinal basis functions on a sequence of hierarchical equidis¬ 
tant nodes are precomputed to machine precision, rendering the interpolation 
problem into a fully parallelisable ensemble of linear combinations of function 
evaluations. A numerical integration algorithm is also proposed, based on inter¬ 
polating the (high-dimensional) integrand. A series of numerical experiments is 
presented, highlighting the practical applicability of the proposed algorithm for 
interpolation and integration for up to 10-dimensional problems. 

In Section[2j we give a brief description of the method, while, in Section[3j we 
show that for the Gaussian kernel the cardinal functions for interpolation can be 
written as a product of univariate cardinal functions. In Section [4j we show that 
if the kernel is a tensor product (such as the Gaussian) then MLSKI is also an 
interplatory scheme. In Section [5j we discuss the integration scheme with some 
numerical results which demonstrate the good approximation properties of the 
method. 

2 Multilevel sparse kernel-based interpolation 

We begin by briefly reviewing the MLSKI method, introduced in m- To con¬ 
struct the sparse kernel-based interpolant, at each level, we solve a number of 
anisotropic radial basis function interpolation problems on appropriate direction- 
wise uniform sub-grids. We briefly note that due to the careful selection of 
anisotropic scalings for the individual interpolation problems, these are typi¬ 
cally sufficiently well conditioned for practical computations; cf. [14] . The indi¬ 
vidually solved interpolation problems are then linearly combined to obtain the 
sparse kernel-based interpolant. The multilevel sparse kernel based interpolation 
(MLSKI) uses a residual interpolation at different levels. 

More specifically, let 17 :=, [0, l] d , d > 2, and let u : 17 —> R. For a multi¬ 
index 1 = £ N d , we define the family of directionally uniform grids 

{Xi : 1 £ N d }, in 17, with meshsize hi = 2 _1 := (2 _il ,..., 2~ ld ). Then, Xi consists 
of the points x^i := (xi u i 17 . .., a n d ,i d ), with = i 3 2~ l \ for i 3 = 0,1,..., 2^, 
j = 1,..., d. The number of nodes TV 1 in Xi is given by 

d 

N 1 = JJ(2*' + 1). 
j=i 

If hi = 2~ n , for all j = 1, • • • , d, Xi is the uniform full grid of level n, having 
size N=(2 n + l) d ; this will be denoted by X n ’ d . 

We also consider the following subset of X” -d+1 ’ d , 

± n ’ d := |J Xi, (1) 

|l|l=n+(d-l) 

with 11|i := Zi + • • • + Id-, which will be referred to as the sparse grid of level n in 
d dimensions ; see Figure^ an illustration with n = 4 and d — 2. 

Further, we define the transformation matrix Ai € by 

A\ :=diag(2 ,1 ,...,2 I <*). 
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Fig. 1: Sparse grid X 4,2 via ([!]). 


The anisotropic radial basis function (RBF) interpolant 5A, of u at the points 
of Xi is then defined by 

Sa,(x) := E ci )i¥ 5(||Ai(x-x M )||), (2) 

Xii GXj 

for x £ 1 ?, where qj £ R. are chosen so that the interpolation conditions 

-SUJx, = «|xi, 

are satisfied. Although (f> can be chosen from a large family of RBF kernels (cf. 
m for details), this work will be focused on the choice of Gaussian kernels, 
that is, cj)(r ) = exp(—c 2 r 2 ), for r > 0 and c > 0 , denoting the, so-called, shape- 
parameter. 

The sparse kernel-based interpolant (SKI, for short) S n ^ is then given by 

E (3) 

9=0 |l|i=n+(<i— 1) — q 

The above combination formula has been used, for instance, for Lagrange poly¬ 
nomial interpolation in [ 8 ], and for the numerical solution of elliptic partial 
differential e qua tions using the finite element method on sparse grids in mm- 
In Figure 2 we show a visualisation of the SKI interpolant S 4 2 as a linear 
combination of the constituent sub-grid interpolants. In the next section, we 
shall show that this procedure indeed results to an interpolation method, when 
4> is Gaussian. 

The convergence of SKI was investigated in m > where it was found that, al¬ 
though it results to acceptable approximation power when interpolating simpler 
tensor-product-type functions /, it can be also prone to very slow convergence 
on more challenging functions. This appears to be due to the fact that, in con¬ 
trast with standard sparse grid/hyperbolic cross/Smolyak-type polynomial or 
linear spline interpolants, the approximation space of the SKI of level n does not 
contain the respective approximation spaces of levels 1 ,..., n — 1 as subspaces. 

In short, the SKI interpolation is not hierarchical in each level. 

To overcome this, a multilevel approach to SKI was proposed in M- There is 
a growing volume of literature on multilevel methods for RBFs, e.g. iisiiflSfliinoisTirflrrai 
. The idea is to interpolate a high frequency residual with a rougher interpolant. 
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Fig. 2: The construction of S 4.2 interpolant on X 4,2 . 


The setting of SKI is naturally suited to be used within a multilevel interpo¬ 
lation algorithm. Firstly, the sparse grids from lower to higher level are nested, 
i.e., X n ’ d C X n+1,d for n £ N. Secondly, each sub-grid interpolant uses appro¬ 
priately scaled anisotropic basis function with the scaling being proportional to 
density of the corresponding constituent sub-grid. Finally, due to the geometri¬ 
cal progression in the problem size from one sparse grid to the next, a multilevel 
algorithm would not affect adversely the attractive complexity properties of SKI. 

The multilevel SKI (MLSKI, for short) algorithm is initialised by computing 
the SKI S n a on the coarsest sparse grid X n o ,d and set A 0 := Sn 0 ,d■ Then, for 
k = 1,.. .n, Ak is the sparse grid interpolant to the residual u — A? on 

X fc,d . The resulting multilevel sparse kernel based interpolant is then given by 

n 

3 =0 

Below we describe how, if the kernel is a tensor product, we can construct 
a tensor product basis for the interpolation once and for all, and so avoid the 
costly solution of subsystems of equations, which is, in principle, highly scalable 
in a parallel computing architecture. 

3 Tensor Products of Univariate Cardinal Functions 

Gaussian kernels on can be viewed as tensor products of univariate Gaussians; 
this is evident also for the anisotropic versions of Gaussians, such as the ones 
discussed in the previous section. This crucial property of Gaussians will be 
instrumental in both the proof of interpolation of the MLSKI algorithm, and the 
development of fast procedures to evaluate the respective MLSKI interpolants. 

We begin by considering the set of cardinal (also known as Lagrange) func¬ 
tions for each interpolation sub-problem in the SKI interpolant Sa,(') hr ([ 3 ]), 
with the particular choice of Gaussian RBFs on the set Xi. Let xi. Xl , be the car¬ 
dinal function on the grid indexed by 1 for the point x^i £ Xi. Let also Xl-,xi. , 
denote the univariate cardinal function in one variable for the node Xi ., with 
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respect to the set of nodes {a?i. . ■ ij = 0,1,, 2 lj }. Hence, for j = 1, ■ • • , d, 
there exist £ M, such that 

2 l J 

(y) = ex P(-c 2 ^( 2 / - 

ij =0 


for y G [0,1]. Then, we have 

d d 

z(y) : = n Xh,x tj , io (yj) = TV ex P(-C 2 h h ( y* - X h,ij ) 2 )> 

j = 1 xuSX] j=i 

where 714 = IIj=i 7 i j} ij ■ This is exactly the form of the cardinal function based 
on the points in the grid Xi, which implies z( y) = \i jXi j, due to the uniqueness 
of Gaussian interpolation. 

Hence, it is possible to compute the cardinal functions for multivariate ap¬ 
proximation by computing a b initio (to arbitrarily high precision, e.g., by using 
symbolic calculators) the cardinal functions for univariate approximation up to 
(for instance) 5,9,17, • • • , 129 equally spaced points, and store these. The ap¬ 
proximation process would then require no solution of linear systems, massively 
increasing the speed of the algorithm. This algorithm will be implemented below 
in the numerical experiments. 

4 Tensor product kernels give interpolatory schemes 

We shall show that the combination formula © is, indeed, an interpolant. To 
highlight the key ideas of the proof, we first consider an example in three dimen¬ 
sions. A two dimensional example is too straightforward, and a four dimensional 
one is already too complicated. 

Setting, for instance, d = 3 and n = 7, we seek to compute the value of 
the sparse kernel-based interpolant 6 ( 3,7 to the function / at the point a = 
(1/4,1/8,1/16) = (2 -2 , 2 -3 , 2 -4 ), which first appears in the grid with multi¬ 
index 1 = (2,3,4). From ([ 3 ]), we see that for d = 3, the sub-grids of the previous 
two levels are linearly combined to give S 37 . Hence, the set of multi-indices 
which index the approximation are {a : |a| = 7, 8, 9}. 

We decompose these sets of multi-indices into two groups: The first sets of 
grids, let us call them Group 1, have two of the three components of the multi¬ 
index less than the corresponding components for the multi-index for the point, 
e.g., the grids represented by the multi-indices {a = (k, 1,1) : k = 5,6,7}. 
The remaining grids, let us call them Group 2, have only one component less 
than the corresponding component for the point; for instance the family of grids 
{a = (k, l, 2) : k + l < 7, k > 2, l > 3}. 

We shall now study the values of the function / on a typical set of grids 
from Group 1, {(/c,l,l) : k = 5,6,7}. To this end, we consider the cardinal 
functions from points on these grids, and their values at the point a. Let b = 
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(ri2 fc , T 2 / 2 , r 3 /2), for some G No, * = 1,2, 3, where 0 < r\ < 2 fc , 0 < r 2 , r 3 < 
2, with cardinal function X(fe,i,i), 6 - Then 

X(fe,l,l),b(a) = Xfc,r 1 2-b(l/4)Xl ) r 2 /2(l/8)Xl 1 r 3 /2(l/16) = 0, 
if ri2 -fc 7 ^ 1/4, since k > 2. If ri2 _fc = 1/4 then 

X(fc,l,l),b( a ) = Xl,r 2 /2(l/8)Xl,r 3 /2(l/16) 
which is independent of k. 

Thus (1/4,r 2 /2,r 3 /2) is a node on all grids {X( fcl>1 ) : k = 5,6,7}, and the 
value of the cardinal functions X(fc,i,i).fc are identical at a. Hence, the contribution 
to the interpolant from /(&), on these grids, is 

2 

EHrO/ww-) 

2 

= /(&)Xl,r 2 / 2 (l/ 8 )Xl,,- 3 / 2 (1/16) ^(-l ) 9 


If we combine the results of the last two paragraphs we see that the contribution 
to the interpolant from all of the points on grids from Group 1 is 0. 

We now turn to grids from Group 2. Let us consider points on the grids X^;^) 
for k + l < 7, k > 2, and l > 3, and their values at a. Let b = (ri2 -fe , r 2 2 _/ , r 3 /4), 
for some 0 < r\ < 2 k , 0 < r 2 < 2 l , and 0 < r 3 < 4 with cardinal function 
X(k,i, 2 ),b- Then, 

X(k,l,2),b(a) = Xfc,r 1 2-fc(l/4)X!,r 2 2-i( 1 /8)X2,r 3 /4(l/16) = 0, 

unless ri2 -fc = 1/4 and T22~ l = 1/8, since k > 2 and l > 3. If ri2 -fe = 1/4 and 
T22~ l = 1/8, then X(k,i, 2 ),b{o) = X2,r 3 /4(1/16), which is independent of k and of 

l. 

Thus, (1/4, l/8,r 3 /4) is contained on all grids X( fei ; i2 ) for k + l < 7, k > 2, 
and l > 3, and the value of the cardinal functions X(k,l, 2 ), 6 > f° r all the permissible 
values of k are identical at a. Hence, the contribution to the interpolant from 
/(&), on these grids, is 

B-D'O") X x<w«,»<«) 

q =0 k-\-l<7—q, k> 2, Z>3 

= /(^)X2,r 3 /4(l/16) y^(-l) g f 2 ) card {(fc, l) : k + l < 7 - q, k > 2, l > 3} 

q= 0 W' 

2 

= /(^)X2,r 3 /4(l/16)^(-l) 9 PV3- g) = 0. 

n \Q/ 




Therefore, the contribution to the interpolant from all grids in Group 2 is also 
zero. 

Thus, the only grid contributing to the interpolant is X( 2 j 3 : 4 ), and the only 
cardinal function from that grid which takes non zero values at a is the cardinal 
function at a itself. Therefore, the SKI S 7,3 is, indeed, an interpolant at all 
points of X 7,3 . Hopefully, this example highlights clearly the key role that tensor- 
product nature of the kernel has in the proof of interpolation property. 

Equipped with the insight gained by the above example, we shall now consider 
the general case. To this end, we begin with the following counting result. 

Lemma 1. Let k G N d and |k| < p G N. Then, 

card {j G N d : j > k, |j| = p} = ^ ^ . 

Proof: For d = 1, the result is immediate, since for k < p, 
card {j G N : j > k, j = p} = 1. 

For d> 2, let us write k = (ki, k), and consider we set , and consider {j G N d_1 : 
|j| = P ~ j 1 }, for ji = fci, ki + 1, • • • ,p- \k\. Then, by induction, 

card {j G N d - X : j > k, |j| = p - ji} = ~ Jl “ d ~ 1 

Thus, 

card {j G N d : j > k, |j| = p} = 


using the well-known summation formula for binomial coefficients; e.g., B 1.2.2]. 

□ 

We are now ready to state and prove the main result of this section. 
Theorem 1. Assuming that the interpolation kernel has the form 

d 

V’(y) = n&ivi), 

i—1 


P~ |k| 

^2 card {j : j > k, |j| = p - ji} 

jl—ki 
P~ |k| 

E 


jl—k 
P-|k|+l 

E 

3 i =i 
p - |k| + d' 
d 


p - ji - |k| + d - 1 

d- 1 


ji+ d—2 

d- 1 


then SKI with this kernel is an interpolatory scheme. 
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Proof: We wish to compute the value of the sparse grid interpolant Sd, n to 
the function / at the point a = (ai2 _Zl , a 2 2 _ * 2 , • • • , ad2~ ld ), with 0 < a 3 - < 2 l \ 
1 < j < d, with at least one of the a 3 £ No odd (this ensures that Xj is the 
grid with the smallest index in which this points appears). From hypothesis, 
the kernel is of tensor-product form, which implies the cardinal functions for 
interpolation are of the form 



as observed in Section [3] 

In order to create a decomposition of the indices as in the example above, 
we need to introduce some notation. Let di,d 2 £ N, with d\ + d 2 = d. Let 
ni G N dl with ni(j) £ {1,2, ••• ,d}, j = 1,2,-■■ ,di, and ni(j) < m(j + 1), 
j = 1,2, • ■ • , d\ — 1. Similarly, let n 2 £ N d2 with n 2 (j) £ {1,2, ••• , d}, j = 
1,2, • • ■ , d 2 , and n 2 (j) < n 2 (j + 1), j = 1,2, • • • ,d 2 - 1. Additionally, rp (j) ^ 
n 2 (fc) for any j, k. In other words, the components of and n 2 exhaust the set 
{1, 2, • • • , d}, and these numbers are ordered within the vectors ni and n 2 . We 
should note that once ni is specified, n 2 is uniquely determined and vice versa. 

Let m, G N di , i = 1,2. Then, let the multi-index m = (ni, mi, n 2 , m 2 ) £ N d 
have components m(fc) = m t (j) if n,(j) = k, i = 1,2. So, for instance, if 
ni(3) = 7, then m(7) = mi(3), and if n 2 (4) = 5, then m(5) = m 2 (4). In this 
way, we break multi-indices into two pieces in a convenient fashion. 

On the other hand, for m £ N d , let m ni £ N dl with m ni (i) = m(ni(i)), 
i = 1,2, ,di and m n2 £ with m„ 2 (i) = m(n 2 (i)), i = 1,2,- ■■ ,d 2 . Then, 
we have the identity 


m = (ni, m ni , n 2 , m„ 2 ). 

Now, for each d\ = 1,2, • • • , d — 1, and ni £ N dl , let 
J(l,di,Hi) = {mj : m! < Z ni }. 


For each mi £ /(l,di,ni), define 

J(l,di,n 1 ,mi) = {(ni,mi,n 2 ,m 2 ) : m 2 £ N d2 , m 2 > l n ,}. 

The sets J(l, d\, ni, mi) partition all multi-indices into sets with a fixed set of 
components of the multi-index less than those of 1 and the remaining components 
greater than or equal to those of 1. The only multi-index missing from this set 
is 1 itself. 

We compute the cardinality of the subsets of J(l, d\, ni, m^, given by 
J (1, di, ni, m 1; g) = {m £ J(l, di,ni, mi) : |m| = n + d - 1 - q}, 
for q = 0, 1, • • ■ , d — 1, which, using LemmajI] is given by 
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We now consider the contribution to the interpolant from a typical point on 
one of the grids indexed by elements of J(l, d\, ni, mi). 

Let m = (n 1 ,m 1 ,n 2 ,m 2 ) € J(l, di, n 1; mi). Then, by definition, l n2 > m 2 . 
Let x = (xi2~ mi , x 2 2- m2 , • • • Xd2~ md ), for some 0 < Xj < 2 mj , i = 1,2, • • • , d, 
with cardinal function Xm.x- Then, 

d 

Xm,x(a) = J__|_ ( a *2 ! ). 

i =1 

Hence, x m ,x( a ) 7^ 0 only if Xi2~ mi = a,2 _Zi , i = n 2 (j), j = 1,2, • • • , d 2 . In this 
case, 


di 

Xm,x(a) = IJ X mn x 2 - m a l( i) («ni(j )2 /ni0) ), 

J=1 

which is independent of m 2 . 

So we have, for fixed ni and mi, the same contribution at a from any indi¬ 
vidual point which appears in a grid in .7(1, d -\, ni, m-|). Thus, the contribution 
to from x at point a is 


9=0 


d — 1 
q 


/(x) Xm,x(a) 

m£J(l,di ,11! ,111! ,q ) 


=/(x) n * 


i=i 




2 m “i(i) 


(a ni (i)2 


^card J(l,d 1 ,n 1 ,m 1 ,g) 
9=0 \ 1 / 


= /(x) fi x. 


i=i 




d -1 


2 -"» 1 ( i) K 1 (i) 2 /ni(5) ) 


xEi- 1 )’ 

9=0 


d-l^^n + d + di -2 — | 1„ 2 1 - q 

d\ — 1 


= 0, 


since the binomial coefficient is a polynomial in <7 of degree less than d— 1 , which 
is thus annihilated by the difference operator. 

Thus, we see that all contributions from points in other grids other than X\ 
are zero at points in X\. Clearly, the only contribution at a G X\ from points in 
X\ will be from a itself. □ 
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5 Numerical Integration 


Equipped with an interpolation method, we can now also discuss integration 
over the unit cube. A quadrature rule will approximate the integral 


[ /(x)dx « V w z /(z), 


where Z C [0, \] d is a finite set of points, and {w z : z £ Z} are the weights. If 
we have a cardinal basis for Z, {\ z : z £ Z}, then we approximate 


f=J 2 f( z ) X: 

zG Z 

giving rise to the quadrature formula 

f(x)dx « /O 


Thus, the weights are given by 


zG z 


[o,i] d 


Xz(x)dx. 


[o.i] d 


Xz(x)dx. 


In the case when the cardinal functions are tensor products, viz., \ z = 
xl, x • • • x Xz d ( as is the case of Gaussians, as observed in Section |3|, we can 
straightforwardly compute the weights as follows: 


p d n 1 

/ Xz(x)dx = TT xL(xj)dXi- 


As also mentioned in Section [3j we can compute the cardinal functions offline 
to arbitrary precision for, e.g., up to 129 equally spaced points, thereby allowing 
for fast and highly accurate computation of the quadrature weights. Indeed, the 
integral of a univariate cardinal function is of the form 


/* 1 r 1 N 

1= x(x)dx = / y; Cj exp(—r 2 (x - yjf)dx, 

Jo Jo i=1 

for some r,yi,Ci £ K, i = 1,- • • , IV, which, upon the change of variables u = 
r(x - 2 /*), gives 


N 




fl ~ r Vi 


i —1 J ~ r ^ 

with the error function erf defined by 


exp (~u 2 )du = ^^2 c *( erf ( r ( 1 “ 2/*) ~ erf (■" r 2/<))> 


N 


2 r x 2 

erf (x) = —= / e~ l dt. 
V 71 " Jo 


The error function can be computed to arbitrary precision. 
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5.1 Examples 


In this section we give a number of examples. We would like to develop a set 
of benchmark examples so we invite other researchers to contact us with their 
results for these cases, and perhaps forward us their own examples. We do two 
tensor product examples, respectively in 5 and 10 dimensions, and then two 
non-tensor product examples, the first being smooth, the second with derivative 
singularities on the boundary. The motivation for the second choice is approxi¬ 
mation of the pay-off function in option pricing problems. 


1. Consider the function 

5 

f{x.) = Y[4xi(l- x i). (5) 

2—1 

The corresponding exact integral of this function on the domain [0, l] 5 with 
16 digits accuracy is (|) 5 = 0.131687242798354. In Table [l] we give the 
MLSKI quadrature computation of the same integral, recording the number 
of nodes use, and the absolute and relative errors. 


nodes 

Absolute error 

Relative error 

243 

3.0091e-2 

2.2850e-l 

1053 

5.1232e-3 

3.8904e-2 

3753 

1.3013e-3 

9.8818e-3 

12033 

1.4927e-4 

1.1335e-3 

36033 

3.6134e-5 

2.7439e-4 

102785 

3.4530e-6 

2.6222e-5 

282525 

8.1811e-7 

6.2125e-6 

754845 

6.9041e-8 

5.2428e-7 


Table 1: MLSKI quadrature for /. 


2. The following example is in 10 dimensions: 

10 

5 (x) = ne^( 1 - )) (6) 

i—1 

The corresponding exact integral of this function on the domain [0, l] 10 with 
16 digits accuracy is 0.194279067580947. Table[2]shows how the error reduces 
with respect to the increase of degrees of freedom. 
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nodes 

Absolute-error 

Relative-error 

59049 

452709 

2421009 

10819089 

1.5068e-l 

5.8153e-3 

3.5882e-3 

4.9348e-4 

7.7556e-l 

2.9933e-2 

1.8469e-2 

2.5400e-3 


Table 2: MLSKI quadrature for g. 


3. Let us now consider the non-tensor product Franke 4D function: 


UF4£>(x) 


iL( —(9^1 — 2) 2 —(9:12 — 2) 2 —(9:13 —2) 2 )/4—(9x 4 — 2) 2 )/8 

4 

i ^ J-(9 xi + 1) 2 )/49-((9x 2 +1) 2 )/10-((9x 3 + 1) 2 )/29-((9x 4 + 1) 2 )/39 

4 

i 1 „(-(9 xi-7) 2 )/4-(9x 2 -3) 2 -((9x 3 -5) 2 )/2-((9x4-5) 2 )/4 

2 

_ 1 e (-(9xi-4) 2 )/4-(9x 2 -7) 2 -((9x 3 -5) 2 )-((9x 4 -5) 2 )_ 


(7) 


The integral of this function on the domain [0, l] 4 with 16 digits accuracy is 
0.037221856819405. The errors in integration are shown in Table [3] 


nodes 

Absolute-error 

Relative-error 

81 

1.6398e-2 

4.4055e-l 

297 

1.2736e-2 

3.4216e-l 

945 

7.9106e-3 

2.1253e-l 

2769 

5.4904e-3 

1.4751e-l 

7681 

5.5825e-4 

1.4998e-2 

20481 

1.3012e-4 

3.4959e-3 

52993 

1.6245e-5 

4.3643e-4 

133889 

1.2027e-7 

3.2312e-6 

331777 

2.2934e-8 

6.1615e-7 


Table 3: MLSKI quadrature for Ufad- 


4. The final example is a 5 dimensional non-tensor product function with 
derivative discontinuities in different directions on the boundary. In order 
to apply MLSKI to derivative pricing problems we would need to be able to 
approximate such functions as 


5 

payoff (x) =^max(a; i 

i— 1 


1 

2 ’ 


0 ). 


( 8 ) 
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The exact integral of this function on the domain [0, l] 5 is |. In Table [ 4 ] we 
see how the error behaves with regard to the degrees of freedom. 


nodes 

Absolute error 

Relative error 

243 

1.5129e-l 

2.4206e-l 

1053 

5.4282e-3 

8.6851e-3 

3753 

2.9705e-3 

4.7529e-3 

12033 

1.0128e-3 

1.6206e-3 

36033 

3.2119e-4 

5.1390e-4 

102785 

9.0693e-5 

1.4511e-4 

282525 

2.2032e-5 

3.5251e-5 

754845 

5.7779e-6 

9.2447e-6 


Table 4: MLSKI quadrature for F payo ff. 
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